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We study the sample-to-sample fluctuations of the overlap probability densities from large-scale 
equilibrium simulations of the three-dimensional Edwards-Anderson spin glass below the critical 
temperature. Ultrametricity, Stochastic Stability and Overlap Equivalence impose constraints on 
the moments of the overlap probability densities that can be tested against numerical data. We 
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I. INTRODUCTION 



Spin glasses are model glassy systems which have been 
studied for decades and have become a paradigm for a 
broad class of scientific applications. They not only pro- 
vide a mathematical model for disordered alloys and their 
striking low-temperature properties (slow dynamics, age- 
dependent response), but they have also been the test- 
ground for new ideas in the study of other complex sys- 
tems, such as structural glasses, colloids, econophysics, 
and combinatorial optimization models. The non-trivial 
phase-space structure of the mean-field solution to spin 
glasses [ll-y] encodes many properties of glassy behavior. 

Whether the predictions of the mean-field solutions 
correctly describe the properties of finite-range spin-glass 
models (and of their experimental counterpart materials) 
is a long-debated question. The Droplet Model describes 
the spin glass phase in terms of a unique state (apart 
from a global inversion symmetry) and predicts a (super- 
universal) coarsening dynamics for the off-equilibrium 
regime. [Jl Moreover, there is no spin glass transition in 
presence of any external magnetic field. On the other 



side, the Replica Symmetry Breaking scenario [3, [5[, 
based on the mean field prediction, describes a complex 
free-energy landscape and a non-trivial order parameter 
distribution in the thermodynamic limit; the dynamics is 
critical at all temperatures in the spin-glass phase. The 
spin glass transition temperature is finite also in presence 
of small magnetic fields; the search for the de Almeida- 
Thoulcss line Tc{h) is the purpose of many numerical 
experiments (see, for example, Ref. 0)- 

From the theoretical perspective, the last decade 
has seen a strong advance in the understanding of 
the properties of the mean-field solution: its correct- 
ness has been rigorously proved thanks to the intro- 
duction of new concepts and tools, like stochastic sta- 
bility or replica and overlap equivalence (THllj. Be- 
sides, numerical simulation has been the methodology of 
choice when investigating finite-range spin glasses, even 
if the computational approach is severely plagued by 
the intrinsic properties (slow convergence to equilibrium, 
slowly growing correlation lengths) of the simulated sys- 
tem's (thermo)dynamics. In this respect, a Moore-law- 
sustained improvement in performance of devices for nu- 
merical computation and new emerging technologies in 



the last years has allowed for very fast-running implemen- 
tation of standard simulation techniques. By means of 
the non-conventional computer Janus [l2| we have been 
able to collect high-quality statistics of equilibrium con- 
figurations of three-dimensional Edwards- Anderson spin 
glasses, well beyond what would have been possible on 
conventional PC clusters. 

Theoretical predictions and Janus numerical data have 
been compared in detail in Refs. [ij and [ij. One of the 
main results presented therein is that equilibrium prop- 
erties at a given finite length scale correspond to out-of- 
equilibrium properties at a given finite time scale. On 
experimentally accessible scales (order 10^ seconds wait- 
ing times corresponding to order 10^ lattice sizes) the 
Replica Symmetry Breaking picture turns out as the only 
relevant effective theory. Theories in which some of the 
fundamental ingredients of the mean-field solutions are 
lacking (overlap equivalence in the TNT model [13, ul- 
trametricity in the Droplet Model) show inconsistencies 
when their predictions are compared to the observed be- 
havior. 

In this work we reconsider the analysis of the huge 
amount of data at our disposal, focusing on the sample- 
to-samplc fluctuations of the distribution of the overlap 
order parameter. The assumptions of the mean-field the- 
ory allow us to make predictions on the joint probabili- 
ties of overlaps among many real replicas which can be 
tested against numerical data for the three-dimensional 
Edwards- Anderson model. The structure of the paper is 
as follows: in section [ll] we give some details on the con- 
sidered spin-glass model and the performed Monte Carlo 
simulations. In the subsequent section we first recall 
some fundamental concepts such as stochastic stability, 
ultrametricity, replica and overlap equivalence and some 
predictions on the joint overlap probability densities, and 
then present a detailed comparison with numerical data. 
In section IIVI we show how finite-size numerical over- 
lap distributions compare to the mean-field prediction 
in which finite-size effects are appropriately introduced. 
We finally present our conclusions in the last section. 



II. MONTE CARLO SIMULATIONS 
A. The Model 

We consider the Edwards- Anderson model [13 in three 
dimensions, with Ising spin variables Ci = ±1 and binary 
random quenched couplings Jy = ±1. Each spin, set 
on the nodes of a cubic lattice of size V = L^ {L being 
the lattice size), interacts only with its nearest neighbors 
under periodic boundary conditions. The Hamiltonian 
is: 



H 



2^ Jtj'^t'^j , 



(1) 



where the sum extends over nearest-neighbor lattice sites. 
In what follows we are dealing mainly with measures of 
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J- mill 


J- max 


Nt 


Ns 


8 


0.150 


1.575 


10 


4000 


16 


0.479 


1.575 


16 


4000 


24 


0.625 


1.600 


28 


4000 


32 


0.703 


1.549 


34 


1000 



TABLE I. A summary of parameters of the simulations we 
have used in this work. For eacli lattice size, L, we considered 
Ns samples, with four independent real replicas per sample. 
For the Parallel Tempering algorithm, Nt temperatures were 
used between Tmin and Tmax, uniformly distributed in that 
range (except in the case of L = 8, in which we have 7 temper- 
atures uniformly distributed between 0.435 and 1.575 plus the 
3 temperatures 0.150, 0.245 and 0.340). Our MCS consisted 
of 10 Heat-Bath sweeps followed by 1 Parallel Tempering up- 
date. More detailed information regarding these simulations 
can be found in Ref. \li 



the spin overlap 






(2) 



where a and b are replica indices, and the sample- 
dependent frequencies Nj{qab) with which we estimate 
the overlap probability distribution Pj{q) of each sample 
(we indicate one-sample quantities by the subscript J): 



Pjiqab) = (siqab-j^J2<''^ 



(3) 



where ((• • ■)) is a thermal average. In what follows (■ ■ ■) 
denotes average over disorder. 



B. Numerical Simulations 

We present an analysis of overlap probability distri- 
butions computed on equilibrium configurations of the 
three-dimensional Edwards-Anderson model defined in 
Eq. (ID). We computed the configurations by means of 
an intensive Monte Carlo simulation on the Janus super- 
computer. Full details of these simulation can be found 
in Ref. Il4l . For easy reference, we summarize the param- 
eters of our simulations in Table U] In order to reach such 
low temperature values, it has been crucial to tailor the 
simulation time, on a sample-by-sample basis, through a 
careful study of the temperature random- walk dynamics 
along the parallel tempering simulation. 



III. REPLICA EQUIVALENCE AND 
ULTRAMETRICITY 

The Sherrington-Kirkpatrick (SK) model T] is the 
mean-field counterpart of model H]). It is defined by the 
Hamiltonian 



H 






(4) 



where the sum now extends to aU pairs of N Ising spins 
and the couphngs Jy- are independent and identicaUy- 
distributed random variables extracted from a Gaus- 
sian or a bimodal distribution with variance 1/A^. The 
quenched average of the thermodynamic potential may 
be performed by rewriting the n-rcplicatcd partition 
function in terms of an nxn overlap matrix Qa.t for which 
the saddle-point approximation gives the self-consistency 
equation 



Qab 



(5) 



where the average {(■ ■ ■)} involves an effective single-site 
Hamiltonian in which Qa,b couples the replicas. The ther- 
modynamics of model (|4]) is recovered in the limit n — >■ 0, 
after averaging over all possible permutations of replicas. 

The overlap probability distribution P{q) is defined in 
terms of such an averaging procedure: for any function 
of the overlap f{q), one has that 



dqa.bP{qa,b)I{qa,b) 



^}^,,-\Y.^^Qpi-)^Pib)) 



(6) 



the sum being over permutations p of the n replica in- 
dices. The assumption of the replica approach is that 
P{q) defined in this way is the same as the large-volume 
limit of the disorder average Pj (q) of the probability dis- 
tribution of the overlap defined in Eqs. ^ and ^. 

The hierarchical solution 0] for Qab is based on two 
main assumptions: stochastic stability and ultrametric- 
ity. In what follows we are interested in the consequences 
of such assumptions when dealing with a generic random 
spin system defined by a Hamiltonian Hj{a), where the 
subscript J summarizes the dependence on a set of ran- 
dom quenched parameters, e.g., the random couplings in 
models ^ and dH). 

Stochastic stability [7|, |8| in the replica formalism is 
equivalent to replica equivalence [3, [l3| : one-replica ob- 
servables retain symmetry under replica permutation 
even when the replica symmetry is broken. This property 
implies that the nxn overlap matrix for an n-replicated 
system, satisfies 



= J2 ifiQac) ~ /(Qbc)] 



(7) 



for any function / and any indices a, 5. In the framework 
of the solution to the mean-field model, this is necessary 
for having a well defined free energy [2, [lO| in the limit 
n ^^ 0. A consequence of (|2l) is, given a set of n real repli- 
cas, the possibility of expressing joint probabilities of m 
among the n(n — l)/2 overlap variables to joint proba- 
bilities for overlaps among a set of up to m replicas. [lO[ 
The following relations hold, for instance, in the cases 



?i = 4, ?Ti = 2 and n = 6, 7ti = 3: 

3P(gi2,g34) = 2P(9i2)P(g34) 

+ ^ (912 - 934) P (912) , (8) 

15P (gi2, 934, 956) = 2P (912, 923, 931 ) 

+ 5P(9)P(9')^(9") 

+ 2<5(g-g')^(9')^(9") 

+ 2<S(g'-9")^(9)P(9') 

+ 2<5(g-g")P(9)P(9') 

+ 2<5(g-g')'5(9'-9")^(9) , (9) 

where q = qi2,q' = 934, 9" = 956- 

Note that relation ([5]) quantifies the fluctuations of the 
overlap distribution: even in the limit of very large vol- 
umes, for the joint probability of two independent over- 
laps. 



P(912,934) = Pj{qi2,q'34) + PjiQi2) PjWzi) ■ (10) 

Ultrametricity is the other remarkable feature of the 
mean-field solution, stating that when picking up three 
equilibrium configurations, either their mutual overlaps 
are all equal or two are equal and smaller than the third. 
A distance can be defined in terms of the overlap so that 
all triangles among states are either equilateral or isosce- 
les. In terms of overlaps probabilities, the property reads: 

^(912, 923, 93l) = (5(912 - 923)'5(923 " 93l)-B(9i2) (11) 

+ [6(912 - 923)^(912, 923)<5(923 - 931) 

-|- two perm.] 

where 0(x) is the Heaviside step function. By replica 
equivalence, A and B can be expressed in terms of 
P{q): M 



A{q,q') = P{q)P{q') , 
B{q) = x{q)P{q) , 

x{q) = r P{q')dq' 



(12) 
(13) 

(14) 



Ultrametricity implies that the joint probability of over- 
laps among n replicas, which in principle depends on 
n{n — l)/2 variables, is a function of only n — 1 vari- 
ables. Thus, using replica equivalence, it is reduced to 
a combination of joint probabilities of a smaller set of 
replicas. Note that P(9i2, 923, 931) is the only non-single- 
overlap quantity appearing in the r.h.s. of Eq. ^: by 
combining replica equivalence and ultrametricity, three- 
overlap probabilities reduce to combinations of single- 
overlap probabilities. 

Stochastic stability, or equivalently replica equivalence, 
is a quite general property that should apply also to 
short-range models, in the hypothesis that the model 
is not unstable upon small random long-range pertur- 
bations [31 • Whether short-range models would feature 
ultrametricity is a long-debated question, for which di- 
rect inspection by numerical means is the methodology 



of choice. It has been shown [23] that, in the hypothe- 
sis that the overlap distribution is non-trivial and fluc- 
tuating in the thermodynamic limit, then ultrametric- 
ity is equivalent to the simpler assumption of overlap 
equivalence, in the sense that it is the unique possibility 
when both replica and overlap equivalence hold. Overlap 
equivalence states that, in the presence of replica sym- 
metry breaking, given any local function Ai{a), the gen- 
eralized overlap qa = iV~^ ^^ ^i(o-")Ai(cr''), with a,b 
indices of real replicas, does not fluctuate when consid- 
ering configurations at fixed spin-overlap [2J|: all def- 
initions of the overlap are equivalent. Assuming that 
stochastic stability is a very generic property, there may 
be violation of ultrametricity only in a situation in which 
also overlap equivalence is violated. In this respect, evi- 
dence of overlap equivalence has been found in both equi- 
librium and off-equilibrium numerical simulations of the 
Edwards- Anderson model [H Hi, [111. 

The aim of this work is a numerical study of the 
samplc-to-samplc fluctuations of the overlap distribution; 
we focus on the sample statistics of the cumulative over- 
lap probability functions defined by 



Xjiq) 



P.J W) dq' 



(15) 




This is a random variable, since it depends on the random 
disorder, and we denote by Ilq{Xj) its probability distri- 



bution. 



We estimate the moments of the 11^ distribution 



Xk{q)= / x''Il,{x)dx=[Xjiq)Y 



Pj{q')dq' 



(16) 



where Pj (q) are the Monte Carlo overlap frequencies for 
a given sample. 

Given a set of three independent spin configurations 
we obtain also the probability for the three overlaps to 
be smaller than q: 



Xr 



J —a 



Pj{qi2,q23,q3i)dqi2dq23dq3i (17) 



In the replica equivalence assumption Xk{q) can be ex- 
pressed in terms of X^i^q) and Xi{q); integrating the 
Ghirlanda-Gucrra relations (|8l9p up to fc = 3 we have: 



(18) 



X,iq) = ^X,{q) + ^XUq) 



Xsiq) 



1 
15 



[2XT{q) + 2Xi{q) 



+ eXfiq) + 5Xf{q)] 



(19) 



Ultrametricity imposes a further constraint: from rela- 
tions (HI]- El) it follows 



XT{q) ^ [xiq)]' = Xtiq) , 



(20) 



FIG. 1. (Color online) The quantity Xt as defined in the text, 
as a function of q for lattice size L = 24 (top) and L — 32 
(bottom) at temperature T ~ 0.64rc. Insets show a magnified 
view of the region q ~ 0.6 (log-log plot). Plots show data for 
Xt computed only with triplets of independent configurations 
(ABC), with triplets in which two configurations belong to 
the same Monte Carlo history (AAB), and triplets in which 
all configurations come from the same Monte Carlo history 
(AAA). No significant difference shows up as long as we take 
enough uncorrelated configurations from the same replica. 



And the quantities ([16)) become polynomials in Xi only. 
The above relation simply states that, if ultrametricity 
holds, the probability of finding three overlaps smaller 
than q factorizes to the probability of finding two overlaps 
independently smaller than q, with the third bound to be 
equal to at least one of the previous two. 

For models in which the overlap is not fluctuating in 
the large- volume limit (i.e., P{q) is a delta function) the 
above relations are satisfied but reduce to trivial identi- 
ties. If the replica symmetry is broken, then stochastic 
stability imposes strong constraints on the form of the 
overlap matrix and consequently on the overlap proba- 
bility densities. Ultrametricity is a further simplification: 
lack of this property might indicate that more than one 
overlap might be needed to describe the equilibrium con- 
figurations [23. 

We can extract further information from the distribu- 



tion Iiq{x). It has been found [26H28J that in mean-field 
theory the probability distribution 7r(y) of the random 
variable Yj ~ 1 — Xj behaves as a power law for Yj ^ 1. 
This implies that Il-q{x) also follows a power law for small 




FIG. 2. (Color online) Top: X2 as a function of the cor- 
responding polynomial in Xi (Eq. H18|l 'l. The straight line 
is the theoretical prediction (unit slope). Center: the ra- 
tio X2/X1 as a function of Xi, where the straight line is 
the theoretical prediction. Bottom: the squared difference 
K2 = [X2 - (Xi + 2X?)/3] ^ as function of Xi . Data refer to 
T ~ 0.64Tc 



X values 



TTg{y ^ I) -- (1 - y) 
n,(s ^ 0) 



x{q)-l 



.;^(«)-l 



(21) 



Since for most samples the Pj{q) is a superposition of 
narrow peaks around sample-dependent q values, sepa- 
rated by wide q intervals in which Pj is exactly zero, 
when dealing with data from simulations of finite-size 
systems, it is convenient to turn to the cumulative dis- 
tribution of the Xj to improve the statistical signal, es- 
pecially at small q values: 



L 


T ~ O.STTc 


T ~ 0.64Tc 


T ~ 0.75rc 


8 


0.625 


— 


0.815 


16 


0.625 


0.698 


0.844 


24 


0.625 


0.697 


0.842 


32 


- 


0.703 


0.831 



TABLE II. Temperature values for each lattice size (Tc = 

1.109 [li,!!^!). 



A. Numerical results 

We recall that in our simulations we tailored the tem- 
perature range for the parallel tempering implementation 
to improve its performance as discussed in Ref. [ij. This 
brought us to direct measurements of observables at tem- 
perature sets that were not perfectly overlapping at all 
lattice sizes. In what follows we compare data at temper- 
atures that are slightly different for different lattice sizes. 
Considering that even if the simulations were performed 
at exactly the same temperatures, tiny size-dependent 
critical effects may always affect the results, we preferred 
not to perform involved interpolations to correct for or- 
der 1% or less of temperature discrepancies. In what 
follows we will refer to the set of data at T ~ 0.64Tc and 
T ^ 0.75Tc for the sake of brevity; the precise values of 
the temperatures are summarized in Table |IT1 We also 
compare data at exactly T = 0.625 = 0.57Tc for lattice 
sizes L = 8, 16,24. 

As our simulations were not optimized to study the 
critical region, we take the value Tc = 1.109(10) from 
Refs. [2^ and [sol (featuring many more samples and small 
sizes to control scaling corrections). Still, combining the 
critical exponents determination of these references with 
the Janus data used herein, we obtain a compatible value 

of 1.105(8). im 

We simulated four independent real replicas per sam- 
ple: thus we avoid any bias in computing XT{q), Eq. (fT7| . 
by picking three configurations in three distinct replicas. 
We show the computed X^iq) for the largest lattices 
L = 24 and L = 32 in Fig. [T] i) considering only con- 
figurations for different replicas (data labeled as ABC); 
ii) picking two configurations out of three from the same 
replica (labeled AAB); iii) picking the three configura- 
tions in the same replica (labeled ^^^) . To minimize 
the effect of bias due to hard samples, we picked up the 
same number of configurations per sample, spaced in time 
by an amount proportional to the exponential autocor- 



n'"(s) = / dxU„{x) (22) relation time r^xp of that sample [Tj]. The three data 

Jo sets (ABC, AAB, AAA) are equivalent and small devia- 



which should verify at small s 



8^(9) 



(23) 



the probability of finding a sample in which the overlap 
probability distribution Pj{q) in the interval [0, q] is small 
enough to verify J^ P{q')dq' < s goes to zero as a power 
law of s. 



sets (ABC, AAB, AAA) are equivalent and small devia- 
tions at low q values remain within error bars: this is a 
strong indication of the statistical quality of our data, as 
described in Ref. [ij. 

We now come to test the Ghirlanda-Guerra relations, 
Eqs. dUl) and (dH). Plotting the two sides of Eq. ^ 
parametrically in q, the data show a slight deviation from 
the theoretical prediction (sec Fig.[2]top). It is interesting 
to compare the discrepancies for different lattice sizes. 
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FIG. 3. (Color online) Data at T ~ 0.64rc. Top: X3 
as a function of the corresponding polynomial in X\ and 
Xt (Eq. mi)). The straight line is the theoretical predic- 
tion (unit slope). Bottom: the squared difference K3 — 
[X3 " (2Xt + 2X1 + 6X1^ + 5Xf )/15]^ as function of Xi, 
T — 0.64rc. Lines connecting points are only a guide to the 
eye. 



FIG. 4. (Color online) Top: The squared difference 

r 91 2 

[Xt— XiJ as a function of Xi. Bottom: the quantity 

K^ = [X3 - {2X1 + SXl + 5X?)/15]^ as a function of X^. 
AU data for T ~ 0.64rc and for lattice sizes L = 16, 24, 32. 
The lines connecting the data points are only intended as a 
guide to the eye. 



As the position and width of P{q) are size-dependent, it 
seems more natural to compare functions of the moments 
Xk for different lattice sizes as functions of the integrated 
probabihty x(g) = Xi{q) (see Fig.[5]middle). It is evident 
from the third plot in Fig. [5] that the quantity 



K2= [X2-{Xi + 2X^)/3y 



(24) 



is definitely non-zero although very small in the entire 
range. However, the data are compatible with K2 de- 
creasing with lattice size and becoming null in the L ^ 00 
limit. 

We can reach similar conclusions regarding X3 as a 
function of X^ and Xi , and the quantity 

K3 = [^3 - (2Xt + 2X1 + 6Xf + 5Xf)/15] ^ (25) 

(see Fig. |3]). Even if the data for different lattice sizes 
stand within a couple of standard deviations, there is a 
clear improvement in the agreement between the predic- 
tion and the Monte Carlo data as the size increases. 

The data plotted in Fig. Intake into account the ultra- 
metric relation ([201) • When comparing Xt and Xf small 
deviations from the prediction arise. However, data for 
L = 32 have strong fluctuations, and do not hint at any 
clear tendency with the system size. The bottom plot in 
Fig. m shows data for the quantity 



A7 = [X3 - {2X1 + 8Xf + 5Xf)/15] ' 



(26) 



which we obtain by substituting pO|) in p5|) . The same 
considerations we made above apply here: the agreement 



with ultrametric relations ((T^ and (|20p improves with 
increasing L. 

Wc can compare the results above with those of Ref . \2&l. 
in which a good agreement between theoretical predic- 
tion of the kind of Eqs. ^, (HH), ^ and Monte Carlo 
data on SD Edwards- Anderson spin glass with Gaussian 
couplings was reported, but without clear evidence on 
whether the very small discrepancies could be controlled 
or not in the limit of large volume. In this respect, we 
have been able to thermalize systems of linear sizes up 
to twice the largest lattice studied in Ref. [2^ and these 
larger sizes show a trend towards satisfying Eqs. (jlSp . 
(fTOll . ([20)1 that was not clear in Ref. [2^. We also note 
that finite-size effects are stronger at low temperatures, 
and obtaining evidence of the correct trend requires data 
from simulations of larger systems than at higher tem- 
perature. We can also compare data at T ~ 0.75Tc and 
T — 0.57Tc (we have data at exactly T = 0.625 for lattice 
sizes L = 8, L = 16 and L = 24 but unfortunately not 
for L = 32). We see that at T - 0.75rc the data for the 
squared differences K^ and [Xt — X^] are almost size- 
independent (this is actually true for [Xt ~ Xf] when 
L > 8, see Fig. El top). At T ~ 0.64rc (see Fig. U), 
such effects cannot be clearly told by comparing only the 
smallest lattices considered, L ~ 16 and L — 24. At 
T = 0.57Tc, size-dependent effects are strong even for 
L = 16, 24 (see Fig. [3 bottom). 

Having data from four independent replicas per sam- 
ple, we have access to the joint probability of two inde- 
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FIG. 5. (Color online) Square difference [Xt - X^] (left) 

and the quantity K^ = [Xg - (2Xi + SXj^ + 5X?)/15]^ 
(right) as a function of Xi. Top: for T = 0.75Tc and 
L = 8, 16, 24, 32. Bottom: for T = O.STTc and L = 8, 16, 24. 



pendent overlaps. According to Eq. ([8]) the quantity 



^(912, 934) 2 D/ I ^ 

-^5(^^^- 3^(912) =P(912|<Z34) 



:^(9 



12 J 



(27) 



(where P{-\-) denotes conditional probability) when plot- 
ted versus qi2, should be a delta function in 534. This 
quantity is shown for L = 32, T ^ 0.64rc and two val- 
ues of 934 in the top plot of Fig. [5] and reveals a clear 
peak around (734. At high 512 values there is a small ex- 
cess in the probability P{qi2) P{<Im), so the difference in 
Eq. ([77]) becomes negative. As one sees in Fig. [5] this 
happens at values gi2 ^ (Zea, i-e., in a region of atypi- 
cally large overlaps that should vanish in the thermody- 
namical limit. The size dependence for the quantity in 
Eq. (P?]) is not easy to quantify from the data: as one 
can see in Fig. [S] (bottom) for a particular choice of 534, 
the peak height tends to increase with L (at least for 
T ^ O.TSTc), but in a very slow way, making extrapola- 
tions in the L — > 00 limit practically impossible. Despite 
this, we note that the negative peaks get narrower as the 
system size increases: we expect then that this effect will 
disappear at larger system sizes. 

We conclude this section commenting the asymptotic 
behavior of the cumulative probability 11^ (z), Eq. (f^. 
The small- 2: decay is clearly a power law (see top plot in 
Fig. [7]), but the best fit exponent is significantly differ- 
ent from the estimate obtained by integrating the overlap 
distribution P{q). Fig. [7] shows a comparison of the ex- 
ponent x{q) obtained by the two methods, for some lat- 
tice sizes, many cut-off values q and two temperatures, 
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FIG. 6. (Color online) Top: The conditioned probability 
^(9121934) (open squares) for L = 32 and T ~ 0.64Tc and 
two values of §34 = 0.211 (left) and 534 = 0.367 (right). We 
also plot 2P(gi2)/3 (open circles) and the difference (full tri- 
angles) of the two above quantities (Eq. H27|) in the text), 
scaled by a factor 2 for a better view. 534 and gsA values are 
indicated by vertical lines for visual reference. We took the 
value <3EA(i = 32, T = 0.64rc) ~ 0.72 as given in Ref. [H. 
Bottom: The difference P{qi2\cisi) — 2P(gi2)/3 with 534 = 
0.367, for different lattice size compared at temperatures 
T = 0.75rc, T = 0.64Tc, T = 0.57rc. 



T - 0.64rc and T - 0.57Tc. Although the differences 
seem to decrease by increasing the lattice size, the trend 
is very slow and even not in a clear direction for some 
values of the cutoff q. Again, the only conclusion that 
can be drawn is that the finite-size effects are large, even 
for L = 32, and safe extrapolations in the L — > 00 limit 
cannot be done. 

A closer inspection of the data reported in Fig. [7] 
reveals that the difference between the two data sets 
is roughly a constant, and this difference becomes ex- 
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FIG. 7. (Color online) Asymptotic behavior of the cumulative 
probability Ug (z) (Eq. ^^). Top: small-z decay for L = 
32, T = 0.64rc and q = 0.3125. Bottom: comparison of 
the exponent x{q) obtained by the two methods described in 
the text (uppermost data points represent values obtained by 
fitting 11,(2), lowermost data points come from integrating 
the P{q)), for some lattice sizes, many cut-off values q and 
temperatures T ~ 0.57rc (left) and T ~ 0.64rc (right). 



L 


T/Tc 


gEA 


a;oo((jEA) 


A 


32 

24 

16 

8 


0.75 
0.64 
0.75 
0.64 
0.57 
0.75 
0.64 
0.57 
0.75 
0.57 


0.663(19) 
0.7319(30) 
0.69674(72) 
0.7625(27) 
0.7954(24) 
0.73780(73) 

0.809(16) 
0.8210(41) 
0.8250(21) 

0.886(18) 


0.91(13) 
0.828(28) 
1.0000(3) 
0.876(24) 
0.842(25) 
1.000031(7) 

1.00(14) 

0.811(49) 

1.000001(9) 

0.95(18) 


0.0923(80) 
0.1015(30) 
0.10618(84) 
0.1182(24) 
0.1216(32) 
0.1443(10) 

0.150(11) 
0.1683(51) 
0.2872(37) 

0.296(28) 


L 


T/Tc 


a 


7 


xVd.o.f. 


32 
24 

16 

8 


0.75 
0.64 
0.75 
0.64 
0.57 
0.75 
0.64 
0.57 
0.75 
0.57 


1.92(34) 
0.93(44) 
2.04(21) 
0.95(21) 
0.75(17) 
1.76(16) 
0.45(21) 
0.53(19) 
0.73(22) 
0.49(16) 


11.2(1.2) 
7.7(1.0) 
9.68(55) 
6.88(41) 
5.62(30) 
5.14(31) 
4.50(52) 
3.37(40) 
2.02(34) 
1.36(17) 


20/97 
38/103 
45/101 
69/107 
88/110 
77/107 
133/113 
161/115 
501/121 
466/123 



TABLE III. Results of the fitting procedure of Eq. ^ on 
numerical P{q) data, with kernel exponent k — 2.5 (see 
Eq. p3|l '). All errors on parameters are jackknife estimates. 

^ We used the symbol x^ in the table to denote the sum of 
squares of residuals, which is not a true chi-square estimator as 
the values of P{q) at different q are mutually correlated. 



IV. THE ORDER PARAMETER 
DISTRIBUTION 



tremely important in the limit of small q, where one 
would expect both measurements of x{q) to approach 
zero. Contrary to expectations, the x{q) estimated from 
the data of 11^ seems to remain non-zero even in the 
(7 — > limit. A possible explanation for this observation 
comes from the fact that the delta peaks in the Pj (q) get 
broader for systems of finite size. Indeed, in the ther- 
modynamic limit, one would expect Pj{q) to be the sum 
of delta functions centered on overlap values extracted 
from the average distribution Pao{q)'- if this expectation 
is true, then the value for Xj{q) is nothing but the prob- 
ability of having a peak at an overlap value smaller than 
q and this is exactly x{q). However, if the delta peaks ac- 
quire a non-zero width A due to finite-size effects, then 
for q < A the overlap probability distribution close to 
the origin Pj{0) may be affected by broad peaks cen- 
tered on overlaps larger than q, which should not count 
in the thermodynamical limit. If this explanation is cor- 
rect, then the limit g -^ for the data shown in Fig. [7] 
(bottom) obtained from n?* should give a rough estimate, 
in the large L limit, for the peak width A (see data in 
Table IIIII and discussion below) . 



We now compare the P{q) obtained in numerical 
simulations of the three-dimensional Edwards- Anderson 
model (dl to the prediction obtained by smoothly intro- 
ducing controlled finite-size effects on a mean-field-like 
distribution consisting in a delta function centered in 
q = qEA and a continuous tail down to q ~ {a simi- 
lar analysis has been carried out for long-range spin-glass 
models, see Ref . 133) . On the positive q axis one has 

Poo(g) = P(g)e(gEA-g) 

+ [1 - a:^oo {qEA)]S{q - ^ea) , (28) 



a;oo(gEA) 



dqPiq) 



(29) 



It is convenient to introduce the effective field h trough 



q ~ tanh (h) 

and consider its distribution 

dq{h) 
'dh' 



(30) 



Voo{h)^Poo{q{h))- 



dqjh) 
dh 



Xoo{q\ 



EA 



^^^ P{q{h))Q{h^p^-h)- 
[1 - Xoo[q-E,A)\5{h - /iea) , 

dhV(h) , 



(31) 
(32) 
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FIG. 8. (Color online) Comparison between the Monte Carlo 
data of the P{q) and the convolution computed as described 
in the text (solid lines). Top: L = 32, T ~ 0.64rc and T ~ 
0.75Tc. Center: L = 24, T ~ 0.57Tc and T ~ 0.64rc. Bottom: 
the conditioned probability P{qi2\q34 ~ go) for L = 32, T ~ 
0.64rc and some values of go- 



being clear that gEA = tanh (/iea)- This change of vari- 
able smooths the constraint on the fluctuations of q near 
the extremes of the distribution. 

In a finite-size system the thcrmodynamical distribu- 
tion Voo{h) will be modified, mainly by the fact that delta 
functions become distributions with non-zero widths. 
Remember that, in the thcrmodynamical limit, we expect 
the distribution 'P.j{h) for any given sample to be the sum 
of delta functions. A simple way to take into account the 
spreading of the delta functions due to finite-size effects 
is to introduce a symmetric convolution kernel 



G^^^(/i-/i') = Cexp -{\h-h'\/A) 



(33) 



where C is a normalizing constant and the spreading pa- 
rameter A is assumed not to depend on /i, Q| while it 
should have a clear dependence on the system size, such 
that limL_j.oo A = 0. The parameter k, to be varied in 



^ This introduces a q-dependent spread, as the Jacobian of the 
transformation (I30I I stretches the distribution at high q values. 



the interval [2, 3], is introduced in order to consider con- 
volutions different from the Gaussian case (fc = 2). 

In order to obtain an analytic expression for the finite 
size distribution 

VL{h) = j dh' ^^^^'^'^^°°^'^'^ Gf{h - h') , (34) 

we assume the following form for the continuous part of 
the distribution 



V{h) ^ P{q{h)) 



dq{h) 
dh 



PiO)il + ah^ + jh^) , (35) 



where ^(O) = -P(O) ~ -Poo(O), a and 7 are free parameters 
to be inferred from the data. The final result is 



VUh) = [I - XooiqEA)] 



Gf{h 



.(k) 



hEA) + G'-^>{h + hEA) 



P(0) 



Ak), 



dz [1 + az' + 72*J G^"^ {h - z) (36) 

''EA 



where Xoo{q-EA) = 2P(0)[/iea + aK^p^/i + lh%Jb]. 

We let a, 7, (/ea and A vary in a fitting procedure to 
P[q) Monte Carlo data; values of P(0) are fixed to the 
Monte Carlo values Pmc{^)- The choice of the exponent 
k in the convolution kernel is crucial. We varied k in the 
interval [2,3]. The Gaussian convolution k — 2 turned 
out to be the worst choice in such interval, giving rise to 
unphysical negative weights for the delta function contri- 
butions, i.e., 1 — a;oo('i'EA) < 0. We obtained very good re- 
sults with the choice k = 2.5. Fit parameters are reported 
in Table Hill for some lattice sizes and temperatures, while 
Fig. |8] shows comparison between Monte Carlo P{q) and 
the relative fitting curve. Although the fitting curves in- 
terpolate nicely the numerical P{q), some of the fitting 
parameters may look strange: in particular gEA is a bit 
larger than the peak location and a;oo((7EA) — 1 (for ex- 
ample, in the L = 32 data the difference is around 2%). 
It is worth remembering that in the solution of the SK 
model at low temperatures the continuous part P{q) has 
a divergence for q — )> g^A' which can easily dominate the 
delta function in finite-size systems (where delta peaks 
are broadened). Indeed, by increasing the system size, 
(7ea seems to move towards the location of the peak max- 
imum and Xoo(?ea) becomes smaller than 1. 

In order to make a stronger test of the above fitting 
procedure, we have used the fit parameters in Table Hill 
to derive the finite-size conditional probability 



PLiq\q')=PL{q,q')/PL{q') 



(37) 



applying the convolution kernel G^ (h~h') to the L — 00 
joint probability given by the Ghirlanda-Guerra relation, 
r.h.s of Eq.ilHI). Fig. [5] shows a comparison between our 
extrapolated PL{qi2\q34 = qo) and the Monte Carlo data 
for L = 32, T = 0.64Tc and three values of qo: the agree- 
ment is very good at any value of qo, especially consid- 
ering that the fitting parameters were previously fixed 
by interpolating the unconditional overlap distribution 
Pdq)- 
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CONCLUSIONS 



We performed a direct inspection of stochastic sta- 
bility and ultrametricity properties on the sample-to- 
sample fluctuations of the overlap probability densities 
obtained by large-scale Monte Carlo simulations of the 
three-dimensional Edwards- Anderson model. Wc found 
small but still sizeable deviations from the prediction of 
the Ghirlanda-Guerra relations but a clear tendency to- 
wards improvement of agreement with increasing system 
size. 

Large fluctuations make it difficult to draw any defini- 
tive conclusion on the analysis of the ultrametric rela- 
tion ([20| when taking into account data for the largest 
lattice size. In addition, critical effects show up at 
T ^ O.TSTc. Considering that for a stochastically stable 
system overlap equivalence is enough to infer ultrametric- 
ity, the results presented here support and integrate the 
analyses and claims of Refs [ij, [ij and [25|, in which the 



authors reported strong evidence of overlap equivalence. 
We also turned our attention to the shape of the 
overlap probability distribution, showing that finite- 
size PL{<i) and PL{q,q') compare well with mean-field 
(infinite-size) predictions, modified by finite-size effects 
that only make delta functions broad. 
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